﻿/*
 * GeographicGeocentricTransform
 * 
 * initial implementation : 30.12.2008
 *
 * Copyright (c) Adam Sida, asida@centrum.cz.
 * All rights reserved.
 * 
 * for license see Petra.Gis.License.txt
 */
using System;
using System.Collections.Generic;
using Petra.Gis.Geometry.Support;


namespace Petra.Gis.ReferenceSystemServices.Transformations
{

  /// <summary>
  /// transform from geographic space (lon, lat) into geocentric (X,Y,Z) forgiven datum
  /// This transform is inversible.
  /// </summary>
  public class GeographicGeocentricTransform : SRSTransform
  {

    // local copy of spheroid parameters variables
    private double a;   // spheroid semi major axis 
    private double e2;  // spheroid fisrt eccentricity
    private double e22; // spheroid second eccentricity


    private double h = 0;

    //Coordinate in3D;
    Coordinate outP;
    Coordinate tmpP;
    AngularUnit angularCsUnit;
    DistanceUnit distanceCsUnit;
    double N, W;




    /// <summary>
    /// xam constr.
    /// </summary>
    public GeographicGeocentricTransform() 
      :base()
    {
      Metadata = new Metadata("Geographic-Geocentric transformation", "none", -1, "", "",
        "transforms geog space lat, lon into 3D (x,y,z) geocentric space");
    }





    public GeographicGeocentricTransform(CoordinateSystem sourceCS, CoordinateSystem targetCS, bool createInverse)
      : base(sourceCS, targetCS)
    {
      IsInverse = createInverse;

      Metadata = new Metadata("Geographic-Geocentric transformation", "none", -1, "", "",
        "transforms geog space lat, lon into 3D (x,y,z) geocentric space");

      UpdateTransform();
    }





    public override bool IsInversible
    {
      get { return true; }
    }






    protected override bool graphCompleted()
    {
      return true;
    }




    protected override void updateTransform()
    {
      //if (_sourceCS == null || _targetCS == null)
      //  return;

      if (!((Source is GeographicCoordinateSystem && Target is GeocentricCoordinateSystem)) ||
            (Source is GeocentricCoordinateSystem && Target is GeographicCoordinateSystem))
        throw new ArgumentException("source and target cs for " + this.GetType().ToString() +
                                    " transform must be one geographic and one geocentric cs");

      GeographicCoordinateSystem tmp;

      if (Source is GeographicCoordinateSystem)
        tmp = (GeographicCoordinateSystem)Source;
      else
        tmp = (GeographicCoordinateSystem)Target;

      angularCsUnit = (AngularUnit)tmp.Unit;
      distanceCsUnit = (DistanceUnit)tmp.Unit;

      a = tmp.Datum.Ellipsoid.SemiMajorAxis * tmp.Datum.Ellipsoid.AxisUnit.MetersPerUnit;
      e2 = tmp.Datum.Ellipsoid.Eccentricity;
      e22 = tmp.Datum.Ellipsoid.SecondEccentricity;

      outP = new Coordinate(0, 0, 0);
      tmpP = new Coordinate(0, 0, 0);

    }





    protected override void transform(ICoordinate P)
    {

      // determine direction
      // if source is geographic and target is geocentric and not inverse
      if (Source is GeographicCoordinateSystem && !IsInverse)
      {
        // convert input lon, lat to geocentric x,y,z
        // convert input lat, lon point to radians
        P.lon *= angularCsUnit.RadiansPerUnit;
        P.lat *= angularCsUnit.RadiansPerUnit;

        W = Math.Sqrt(1 - e2 * Math.Pow(Math.Sin(P.lon), 2));
        N = a / W;

        // lon,lat to x1,y1,z1
        outP.X = (N + h) * Math.Cos(P.lon) * Math.Cos(P.lat);
        outP.Y = (N + h) * Math.Cos(P.lon) * Math.Sin(P.lat);

        if (P.Z == 0)
          outP.Z = (N * (1 - e2) + h) * Math.Sin(P.lon);
        else
          outP.Z = (N * (1 - e2) + P.Z) * Math.Sin(P.lon);

        //convert to output cs units
        P.X = outP.X / distanceCsUnit.MetersPerUnit;
        P.Y = outP.Y / distanceCsUnit.MetersPerUnit;
        P.Z = outP.Z / distanceCsUnit.MetersPerUnit;
      }
      else
      {
        // inverse = convert geocentric x,y,z to lon, lat
        // convert input x, y, z do meters
        P.X *= distanceCsUnit.MetersPerUnit;
        P.Y *= distanceCsUnit.MetersPerUnit;
        P.Z *= distanceCsUnit.MetersPerUnit;

        // x2,y2,z2 to lon,lat
        outP.lat = Math.Atan(P.Y / P.X);
        outP.lon = Math.Atan(P.Z * (1 + e22) / Math.Sqrt((P.X * P.X) + (P.Y * P.Y)));

        W = Math.Sqrt(1 - e2 * Math.Pow(Math.Sin(outP.lon), 2));
        N = a / W;

        // iteration to get highest possible precision
        for (int t = 0; t <= 5; t++)
        {
          outP.lon = Math.Atan((P.Z + N * e2 * Math.Sin(outP.lon)) / Math.Sqrt((P.X * P.X) + (P.Y * P.Y)));
          W = Math.Sqrt(1 - e2 * Math.Pow(Math.Sin(outP.lon), 2));
          N = a / W;
        }

        // convert to output angular unit
        P.lon = outP.lon / angularCsUnit.RadiansPerUnit;
        P.lat = outP.lat / angularCsUnit.RadiansPerUnit;

        // prepare height...
        // todo: height must be converted to output coordinate system vertical datum units
        // if no vertical datum defined z remains in meters as now
        W = Math.Sqrt(1 - e2 * Math.Pow(Math.Sin(outP.lon), 2));
        N = a / W;

        tmpP.X = (N + h) * Math.Cos(outP.lon) * Math.Cos(outP.lat);
        tmpP.Y = (N + h) * Math.Cos(outP.lon) * Math.Sin(outP.lat);
        tmpP.Z = (N * (1 - e2)) * Math.Sin(outP.lon);

        outP.Z = Math.Sqrt(Math.Pow(P.X - tmpP.X, 2) + Math.Pow(P.Y - tmpP.Y, 2) + Math.Pow(P.Z - tmpP.Z, 2));

        P.Z = outP.Z;
      }
    }






    public override SRSTransform Inverse()
    {
      return new GeographicGeocentricTransform(Source, Target, !IsInverse);
    }





    public override Dimension TransformDimension
    {
      get { return Dimension.dim3D; }
    }


  }
}
